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Abstract. Observations of submillimeter lines of CO, HCO + , HCN and their isotopes from circumstellar disks 
around low mass pre-main sequence stars are presented. CO lines up to J=6— >5, and HCO + and HCN lines up 
to J=4^3, are detected from the disks around LkCa 15 and TW Hya. These lines originate from levels with 
higher excitation temperatures and critical densities than studied before. Combined with interferometer data on 
lower excitation lines, the line ratios can be used to constrain the physical structure of the disk. The different 
line ratios and optical depths indicate that most of the observed line emission arises from an intermediate disk 
layer with high densities of 10 6 — 10 s cm -3 and moderately warm temperatures in the outer regions. The data are 
compared with three different disk models from the literature using a full 2D Monte Carlo radiative transfer code. 
The abundances of the molecules are constrained from the more optically thin 13 C species and indicate depletions 
of ~ 1 — 30 for LkCa 15 and very high depletions of > 100 for TW Hya with respect to dark cloud abundances. 
Evidence for significant freeze-out (factors of 10 or larger) of CO and HCO + onto grain surfaces at temperatures 
below 22 K is found, but the abundances of these molecules must also be low in the warmer upper layer, most 
likely as a result of photodissociation. A warm upper layer near the surface of a flaring disk heated by stellar and 
interstellar radiation is an appropriate description of the observations of TW Hya. LkCa 15 seems to be cooler at 
the surface, perhaps due to dust settling. The density constraints are also well fitted by the flared disk models. 
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1. Introduction 

Circumstellar disks play an essential role in the under- 
standing of the formation of planetary systems such as 
our own (see Beckwith 1999 and Beckwith & Sargent 
1996 for recent reviews). These protoplanetary disks con- 
tain a few percent of the mass of the pre-main se- 
quence stars which they surround. One of the key ques- 
tions concerning circumstellar disks is their evolution to- 
ward planetary formation. The different evolution scenar- 
ios can be constrained by placing limits on the density 
and temperature distributions in the disks. The standard 
method for determining the disk physical structure utilizes 
fits to the observed spectral energy distributions (SEDs) 
( Adams, Shu fc Lada 1988| ). This procedure relies on the 
changing opacity of the dust at different wavelengths. At 
long wavelengths (typically A >1 mm), the dust emis- 
sion is optically thin and hence traces the product of 
mass and mean temperature (Beckwith 1999), whereas 
at shorter wavelengths the disk becomes optically thick 
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so that only the temperature structure and geometry of 
the disk surface-layer is probed. The derived temperature 
and density solution is not unique since different distribu- 
t ions or different dust properties are able to fit the SEDs 
( [Bouvier fe Bcrtout 19"9^ , |Thamm et al. 1994| . In addi- 
tion, changing dust properties with position in the disk 
can affect the analysis, as can the disk size. Nevertheless, 
one of the more robust results has been the recognition 
of relatively high temperatures in the surface layers of the 
disk, implying that they need to be heated more efficiently 
by stellar radiation compared to the traditional thin (flat) 
disk model. This led Kcnyon fc Hartmann (1987)] to pro- 
pose a flared disk geometry in which the outer disk inter- 
cepts more radiation than does a flat disk. 

Hubble Space Telescope (HST) observations of young 
low mass stars such as HH 30 and HK Tauri show 
edge on (silhouette) disks which indeed flare noticeably 
( Burrows ct al. 1996 ). The radiation from the central star 
incident on the outer parts of the disk changes the tem- 
perature and chemistry in those regions, with the tem- 
perature change giving rise to a larger scale height and 
thereby flaring the disk. Recent models by Chiang & 
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Goldreich (1997, 1999) and D'Alessio et al. (1997, 1998, 

1999) include the irradiation of flared disks to derive self- 
consistent models with a warm outer layer. The models 
by Bell et al. (1997, 1999) take both the stellar radiation 
and re-processing of radiation in the disk into account. 
The latter models have an isothermal temperature in the 
vertical z— direction due to large flaring in the inner disk 
region, thereby shielding the the outer disk from stellar 
light. Comparison with the other models provides a good 
test case whether a high temperature upper layer is needed 
to satisfy the observational constraints. All three types of 
models are used in this work and will be discussed in more 
detail in §|. 

An alternative method to derive the density and tem- 
perature structure in disks is through modeling of molec- 
ular line emission. Although the inferred solution from 
observations of a single line is not unique, data on a suf- 
ficiently large number of transitions of various molecules 
can be used to constrain the temperature and density in- 
dependently. Moreover, careful analysis of the line profiles 
can provide positional information, since the center of a 
line probes a different radial part of the disk compared 
with the wings, unless the disk is nearly face-on. In addi- 
tion, observations of various isotopomers can give informa- 
tion on different vertical regions of the disks due to their 
varying optical depths. To date, most data concern the 
lowest rotational J=l-0 and 2-1 transitions of 12 CO and 
13 CO, which originate from levels at low energies (< 20 K) 
and which have low critical densities (< 5000 cm -3 ) (e.g., 
Dutrey et al. 1996). Data on molecules with larger dipole 
moments such as HCN and HCO + have been limited to 
the 1.3 millimeter band (Dutrey et al. 1997), except for the 
case of TW Hya (Kastner et al. 1997). In this paper, higher 
rotational lines in the 0.8 and 0.45 millimeter atmospheric 
windows are presented, obtained with the James Clerk 
Maxwell Telescope (JCMT) and Caltech Submillimeter 
Observatory (CSO). These lines probe higher tempera- 
tures (up to 100 K) and higher densities (up to 10 7 cm~ 3 ) 
than do presently available spectra. 

The observations are accompanied by a detailed anal- 
ysis of the excitation and radiative transfer of the lines. 
In contrast with previous models (e.g. Gomez & D'Alessio 

2000) , our analysis uses statistical equilibrium (SE) rather 
than local thermodynamic equilibrium (LTE) since the 
surface layers of the disk may have densities below the 
critical density of various transitions. In addition, the two- 
dimensional (2D) code developed by Hogerheijde & van 
der Tak (2000) is used to calculate the full radiative trans- 
fer in the lines. The data can be used to test the disk 
models described above that are fit to the SEDs available 
for most T-Tauri and Herbig Ae stars. In addition to con- 
straining the temperature and density, the observations 
and models also provide information on the depletion of 
different species. 

The molecular abundances and excitation are studied 
by comparing different isotopomers of CO, HCO + and 
HCN for two sources: TW Hya and LkCa 15. TW Hya 
is nearby (57 pc, Kastner et al. 1997) and has a disk 



seen nearly face-on. LkCa 15 is located at the edge of the 
Taurus cloud at ^140 pc and has an inclination of ~60°, 
where 0° is face-on. Both sources show a wealth of molec- 
ular lines and are well-suited for developing the analysis 
tools needed to investigate disk structure. 

The outline of this paper is as follows. In §2, we present 
the observational data. In §3, we perform a simple zeroth- 
order analysis of the observed line ratios to constrain the 
excitation parameters. The adopted disk models are in- 
troduced in §4.1, whereas the methods for calculating the 
level populations are explained in §4.2-4.5. Finally, the 
results of the analysis are given in §5 and summarized in 
§6. 

2. Observations 

Between September 1998 and December 2000, spec- 
tral line observations were obtained for several pre- 
main sequence low mass stars surrounded by circum- 
stellar disks using the dual polarization B3 receiver at 
the James Clerk Maxwell Telescope (JCMT^|) in the 
345 GHz (0.8 mm) band. The observations were ob- 
tained mostly in single side band (SSB) mode us- 
ing beam-switching with a typical switch of 180" in 
azimuth. The spectra were recorded with the Digital 
Autocorrelation Spectrometer (DAS) at a frequency reso- 
lution of ^0.15 MHz (~ 0.15 km s _1 ), and were converted 
to the main-beam temperature scale using ?7 m b=0.62. See 



< ittp : //www . j ach . hawaii . edu/ JACpublic/ JCMT/| > for 

details. The calibration was checked regularly at each fre- 
quency setting against standard spectra of bright sources 
obtained by the JCMT staff, and were generally found to 
agree within 10%. Integration times ranged from 30 min- 
utes (ON+OFF) for 12 CO 3-2 up to 120 minutes for C 18 
3-2, reaching rms noise levels on the T m t, scale of about 
20 mK after adding the data from the two mixers together 
and smoothing to 0.3 MHz resolution. A deep integration 
on the C 18 2-1 line was obtained with receiver A3 for 
LkCa 15, which has ry m b = 0.79. 

These data are complemented by observations using 
the Caltech Submillimeter Observatory (CSO^J) of the 
12 CO 6-5 line for the same sources. In addition, the 12 CO 
2-1 line has been observed with the IRAM 30m telescope^ 
for LkCa 15. For the 12 CO 6-5 line, ?7 m b=0.40, whereas 
for the IRAM 12 CO 2-1 observations the raw data are 



divided by 0.55 (see < ittp : //www . iram. es/ >). 

Interferometer maps of the lowest rotational transi- 
tions of several species toward LkCa 15 have been ob- 



tained by Qi (2000) using the Owens Valley Millimeter 
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Fig. 1. Top: selected CO, HCO + and HCN lines toward LkCa 15. The profiles show a double-peaked structure typical for a 
disk seen at an inclination of about 60°. Bottom: selected CO, HCO + and HCN lines from the face-on disk around TW Hya. 



Array (OVRO). Some lines have also been imaged by 
Simon et al. (2000) and Duvert et al. (2000) with the 
IRAM Plateau de Bure interferometer. In addition, the 
Infrared Space Observatory (ISO) has detected the lowest 
rotational S(0) and S(l) lines of H 2 , which provide inde- 
pendent constraints on the temperature and mass of warm 
gas and which are discussed elsewhere ( Thi ct al. 2001 ) . In 
this paper only the single dish results on CO, HCO + and 
HCN for the sources LkCa 15 and TW Hya are presented. 

Figure [l] shows some of the spectra observed toward 
LkCa 15 and TW Hya. The double peaked profiles for 
LkCa 15 are consistent with Keplerian rotation of the 
disk seen at an inclination of 58 ± 10° (Qi 2000, Duvert 
et al. 2000). Since TW Hya is seen face-on, only narrow 
single-peaked lines are observed from this source. For both 
stars, the 12 CO lines disappear at one beam offset from the 
source. Table 1 summarizes the measured line parameters 
and beam sizes at the observed frequencies. The upper- 
limits for LkCa 15 refer to a 2xrms noise level, with the 
limit on the integrated line strengths obtained by using 
two separate gaussians each with a line- width of 1.3 km 
s" 1 , as found for 13 CO 3-2. For TW Hya, the upper limits 
assume a gaussian with a width of 0.76 km s _1 , similar to 



that observed for HCN and H 13 CO+ 4-3. Note that our 
HCO + 4-3 line toward TW Hya is a factor of three weaker 
than that found by Kastner et al. (1997). We adopt our 
values in the analysis. The HCN 4-3 integrated intensity is 
comparable to that found by Kastner et al. (1997) within 
10 %. There is a hint of a 12 CO 6-5 line toward TW Hya, 
but this is treated as an upper limit. 



3. Simple one-dimensional analysis 

Although the observed line intensities are a complex 
function of the physical structure of the disk and the 
line/continuum optical depth, useful insights can be ob- 
tained from a simple one dimensional analysis of the line 
ratios. For constant temperature and density models such 
as presented by Jansen et al. (1994) and Jansen (1995), 
the data provide constraints on both parameters. To com- 
pare data obtained with different beams the intensities 
were scaled to the same beam (see §4.3). 

Consider first the observed ratios of 12 CO and its iso- 
topomers. The 12 CO 3-2/ 13 CO 3-2 ratios of 3.3 and 7.6 
for LkCa 15 and TW Hya, respectively, indicate that the 
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Table 1. Observed line parameters for LkCa 15 and TW Hya 
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j T mb dv 
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Date 




K km s" 1 
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CO 6-5 


0.53 


0.29/0.28 
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14.5 


CSO 


Jun '00 


CO 3-2 


1.39 


0.60/0.56 


3.3 


13.8 


JCMT 


Nov '99 


CO 3-2 


1.17 


0.37/0.39 


2.2 


25.7 


CSO 


Feb '98 


CO 2-1 


1.82 


0.74/0.76 


2.9 


10.5 


IRAM 


Dec '98 


13 CO 3-2 


0.39 


0.13/0.15 


3.4 


14.4 


JCMT 


Sept '98 


13 CO l-0 e 


7.43 






3.1x2.6 


OVRO 




C ls O 3-2 


<0.14 6 


<0.05 c 




14.5 


JCMT 


Nov '99 


C 18 2-1 


<0.20 b 


<0.07 c 




22.2 


JCMT 


Nov '98 


C 18 l-0 e 


0.70 






4.3x4.0 


OVRO 




HCO+ 4-3 


0.26 


0.14/0.14 


3.3 


13.4 


JCMT 


Sept '98 


HCO+ l-0 e 


4.19 






4.5x3.3 


OVRO 




H 13 CO+ 4-3 


<0.13 b 


<0.05 c 




13.7 


JCMT 


Jan '00 


H 13 CO+ l-0 e 


0.07 






13.0x10.8 


OVRO 




HCN 4-3 
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0.09/0.08 


3.3 


13.5 


JCMT 


Sept '98 


HCN l a - Of 


3.04 






4.3x3.4 


OVRO 




H 13 CN 3-2 e 


1.49 






0.9x0.6 


OVRO 




H 13 CN 1 2 - Of 


1.20 






5.8x4.6 


OVRO 




TW Hya 


CO 6-5 


< 3.22 


< 1.19 


2.46 


14.5 


CSO 


Jun '00 


CO 4-3 d 


5.0 






11.0 


JCMT 




CO 3-2 


1.98 


2.94 


0.63 


13.8 


JCMT 


Nov '99 


CO 3-2 


1.00 


0.77 


1.23 


25.7 


CSO 


Jun '00 


CO 2-l d 


1.02 






20.0 


JCMT 




13 CO 3-2 


0.24 


0.29 


0.78 


14.4 


JCMT 


Feb '99 
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0.14 
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13.5 
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HCN 3-2 d 
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a Width of best single gaussian fit to total profile 

b Calculated assuming the line is double peaked consisting of two separate gaussians, each with a width of 1.3 km s~ 



Listed value is 2xrms 

Values from Kastner e t al. (1997) 



e Values from |Qi (2000) 

f Calculated assuming a line width of 0.76 km s" 1 



CO lines are optically thick, assuming a normal isotope 
ratio of [ 12 C]/[ 13 C]=60 in the solar neighborhood. On the 
other hand, the 13 CO emission has an optical depth of 
only a few, since the C 18 3-2 emission is not detected. 
The ratio of peak 3-2 temperatures of >3 (using the 2a 
limit for C 18 0) and the observed beam-corrected 1-0 ratio 
of 5.0 are only marginally smaller than the isotope ratio 
[ 13 C]/[ 18 0] of 8.3 (Wilson & Rood 1994). 

The temperature can be determined from the 13 CO 
3-2/1-0 ratio of 1.35±0.4, which gives temperatures of 
~20-40 K in LkCa 15. Care should be taken with the in- 
terpretation of this result since the emission of the two 
lines most likely comes from different regions of the disk 
due to the difference in optical depth of the two lines (see 
§|3). The beam-corrected ratio for 13 CO 3-2/2-1 of 0.9 
for TW Hya indicates that the bulk of the gas in this 
source is colder than 25 K for densities > 10 5 cm -3 . The 



12 CO 6-5 line probes higher temperatures since its upper 
level lies at an energy of 116 K. The observed ratio 12 CO 
6-5/3-2=0. 42l°' 2 4 for LkCa 15 also suggests the presence 
of gas with a temperature in the range 20-40 K while 
the upper limit of CO 6-5/ 3-2 < 1.0 for TW Hya gives 
T <150 K (cf. Figure 4 of |Janscn et al.l996j ). The 4-3 
(JCMT) / 3-2 (CSO) ratio of 0.91+^ suggests temper- 
atures ~40 K whereas the ratio of both JCMT lines indi- 
cates somewhat higher temperatures. The different optical 
depths of the 3-2 and 6-5 lines imply that they probe dif- 
ferent vertical layers of the disk. Such vertical structure 
may affect these conclusions, though not by large factors 
(see fO). 



The density is best probed by molecules with a large 
dipole moment such as HCO + and HCN. The measured 
HCO+/H 13 CO+ 1-0 ratio of 6.2 for LkCa 15, and the 4-3 
ratio of 6.7 for TW Hya and >1.9 for LkCa 15 indicate 
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that the main isotopomeric lines are again optically thick. 
No limits on HC 18 0+ exist, but H 13 CO+ may be close to 
optically thin (see §jO). The limit on the H 13 CO+ 4-3 



line toward LkCa 15 together with the 1-0 line detected 
with OVRO gives a 4-3/1-0 ratio of less than < 2.4 and 
constrains the density to < 10 7 cm -3 in the HCO + emit- 
ting region. The optically thick HCO + 4-3/1-0 ratio of 

0. 751.;i| suggests n > 10 6 cm" 3 at T =20-30 K. The 
HCN l 2 -0i to H 13 CN l 2 -0i ratio of 1.4 indicates that 
both lines are severely optically thick. The HCN 4-3 line 
has an even higher critical density than that of HCO + 4- 

3. The observed HCN 4-3/1-0 ratio of l-0±g | indicates 
densities n « 10 7 — 10 8 cm~ 3 in the LkCa 15 disk. 

For the TW Hya disk, the HCN/H 13 CN 4-3 ratio has 
a lower limit of 12.3, indicating that the HCN lines are 
optically thin or nearly optically thin. The HCN 4-3/3-2 
ratio of O.BIq' 3 , constrains its density to lie in the 10 6 to 
10 s cm -3 range, and, as the lines are (nearly) optically 
thin, this should refer to the regions in the disk where 
HCN is most abundant. 

In summary, the simple analysis indicates that the 
main isotope lines are optically thick, but that the lines of 
13 C isotopomers of CO and HCO + have at most moderate 
optical depths. The bulk of the gas is cold, but the pres- 
ence of a warm layer is suggested from the CO 4-3 line for 
TW Hya. The inferred densities in the region where the 
lines originate are high, at least 10 6 cm -3 , but not suf- 
ficient to thermalize all transitions, especially those from 
high dipole moment species. 

4. Description of models 

4.1. Adopted disk models 

In this work, the line emission from three recently pub- 
lished disk models is calculated and compared with ob- 
servations. Although each of these models has limitations, 
they are representative of the range of temperatures and 
densities that may occur in disk models, even if not specifi- 
cally developed for the large radii probed in this work. The 
three disk treatments analyzed in detail arc: 

1. The model by D'Alessio et al. (1999) (see also 
D'Alessio, Calvet & Hartmann 1997 and D'Alessio et 
al. 1998), in which the disks are geometrically thin and 
in vertical hydrostatic equilibrium. The gas and dust 
are heated by viscous dissipation, radioactive decay, 
cosmic rays and stellar irradiation. The gas and dust 
are coupled and can thus be described by a single tem- 
perature. The assumption of a geometrically thin disk 
implies that the energy balance can be calculated at 
each radius separately, decoupled from other radii. The 
temperature and density distribution (Fig. ^, top two 
Figs.) have been calculated using this procedure. 

2. The model by Chiang & Goldreich (1997, 1999) repre- 
sents a passive disk in radiative and hydrostatic bal- 
ance. The disk structure is bi-layered with a super- 
heated dust surface layer that is in radiative balance 
with the light from the central star. This super-heated 
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Fig. 2. Comparison of the density in cm -3 (left three figures) 
and the temperature in K (right three figures) for the models 
of D'Alessio et al. (1999) (top), Chiang & Goldreich (1997) 
(middle) and Bell (1999) (bottom). All models have a gas+dust 
mass of 0.024 Mq. Only one quadrant of the 2-D flaring disk 
is shown. 



layer radiates half of its energy toward the midplane, 
thereby heating the inner regions of the disk. Figure 
^| (center two figures) presents the temperature and 
density distribution from this model. 
3. The model by Bell et al. (1997) and Bell (1999) takes 
the effect of both stellar irradiation and the reprocess- 
ing of radiation into account. This redistribution of en- 
ergy gives rise to a strong flaring of the disk near the 
star, but less strongly further away. The assumption of 
vertical hydrostatic equilibrium is used; however, the 
flaring of the disk at small radii (the first few AU) may, 
for high values of the mass accretion rate M, shield the 
surface of the disk at larger radii from the radiation. 
The resulting temperature and density distribution is 
shown in Figure || (bottom two figures). The shield- 
ing of the outer regions results in a cold disk which is 
approximately isothermal in the vertical direction. 

The published models described above were calculated 
for different disk masses. Since the results depend sensi- 
tively on mass and cannot simply be scaled, this greatly 
complicates their comparison. To eliminate this uncer- 
tainty, the authors have kindly supplied several of their 
model runs for various masses. The model presented here 
for LkCa 15 has a mass of 0.024 Mq, close to its ob- 
served value (Osterloh fc Beckwith 1995). The m 
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the TW Hya disk estimated from submillimeter and cen- 
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timeter continuum observations is ^0.03 M Q , assuming a 
gas/dust ratio of 100/1 (Holland ct al. 2000, Wilner et 
al. 2000). The disks all have an outer radius of 200 AU, 
and the precise mass is fixed by dividing the density of 
an appropriate model by a small factor of less than 3. 
Figure [5] compares the fraction of mass at a given tem- 
perature or density in the three models, whereas the dis- 
tributions within the disk are shown in Figure |2|. While 
the density distributions are similar, the temperature dis- 
tributions between the three models are clearly different. 
One of the aims of this paper is to investigate whether 
the molecular line data are consistent with these different 
types of models. 

The adopted models are not tailor-made for the two 
sources studied here. For example, they do not fit in de- 
tail the observed SEDs: the Bell (1999) model is too cold 
on the outside to reproduce the mid-infrared emission. 
Chiang et al. (2001) have presented models for LkCa 15 
and TW Hya which fit the observed SEDs. However, both 
models have significant setting of the dust. The gas may 
still flare out to higher vertical distances, but the SED 
does not provide observational constraints. For this rea- 
son, we adopted the original Chiang & Goldreich (1997) 
model which has no settling of dust so that the temper- 
ature is defined over the entire disk. Other parameters 
entering the models are the disk accretion rate and the 
luminosity and effective temperature of the star. For the 
accretion rate, which enters the D'Alessio et al. models, a 
value of 10" 8 M yr" 1 and a = 0.01 was chosen, which 
is higher than the observed values of 10 -9 and 5 x 10~ 10 
M Q yr _1 fo r LkCa 15 and TW Hy a, respectively (Hartman 
et al. 1998, Muzcrollc ct al. 2000| ). However, the observed 
molecular lines probe the outer region of the disk whereas 
the accretional heating due to the high values of a and 
the accretion rate will only affect the inner few AU. The 
models refer to a 0.5 M Q star with T cff =4000 K. For 
comparison, LkCa 15 has a mass of roughly 1 M Q and 
an effective temperature of 4400 K (Siess et al. 1999), 
while TW Hya has a mass of 0.7 M and T eff = 4000 K 



(Muzerolle et al. 200C) 



4.2. Radiative transfer methods 

The radiative transfer in the molecular lines from disks 
is calculated in two steps. First, the abundances of the 
molecules in the disk are estimated using a ray-tracing 
method in the vertical direction through the disk and 
adopting statistical equilibrium calculations. The ratio of 
the different modeled lines constrains the range of deple- 
tions. This calculation does not take into account the incli- 
nation of the source and assumes that the ratio of different 
lines is less sensitive to inclination than the integrated in- 
tensity of a single line. 

Once the depletions are constrained, a full 2D radia- 
tive transfer calculation is performed using the acceler- 
ated Monte Carlo (AMC) code of Hogerheijde & van der 
Tak (2000), whose results are compared to the observa- 
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Fig. 3. Comparison of the fraction of mass in a given temper- 
ature and density interval for the models of D'Alessio et al. 
(1999) (top figures), Chiang & Goldreich (1997) (middle fig- 
ures) and Bell (1999) (lower figures). The distributions in the 
disks are plotted in Figure [| 

tions. The motivation for this elaborate scheme is given 
in §4.3 and is driven by the huge computational time in- 
volved in the latter calculation. The AMC code has been 
compared with other radiative line transfer codes in a 
workshop in Leiden in 1999, where the populations of 
the levels and convergence have been tested for a set of 
one-dimensional problems. The comparison of the various 
codes is described in van Zadelhoff et al. (in preparation). 



4.3. Level populations and depletions 

The level populations x u and xi can be calculated in var- 
ious ways. The simplest approximation is that of Local 
Thermodynamic Equilibrium (LTE), which is valid for all 
levels that are collisionally excited in a gas with densities 
higher than the critical density for that level. The latter 
is given by 



A, 



(1) 



where J2i Kui denotes the sum of all collisional rate co- 
efficients from level u to all other levels i. LTE is usually 
adopted in the analysis of line emission from disks (e.g., 
Gomez & D'Alessio 2000, Dutrey et al. 1996, Kastner et 
al. 1997). Although the bulk of the gas is in LTE due to the 
high densities, the radiation from the central part of the 



see also http: / /www.strw.lcidenuniv.nl/^radtrara 
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Fig. 4. The relative populations of levels of the CO molecule 
(J=3 and J— 6) and the HCO + molecule (J=l and J— A) for 
three different calculations (LTE (dotted), SE with no stim- 
ulated emission or absorption (dashed) and SE with full 2D 
radiative transfer (solid)). The levels are calculated for the 
D'Alessio et al. (1999) model and plotted as a function of height 
Z at a radius of 175 AU. The assumed abundances are 10 -4 
for CO and 5-10" 9 for HCO+. In the lower plot the HCO+ 
J=4 level population is reduced by a factor of 5 for clarity. It 
is clearly seen that the populations in the upper layer of the 
disk are not in LTE. 

disk can be highly optically thick depending on the molec- 
ular abundances and assumed depletions. In that case the 
outer layers dominate the emission where the densities 
can fall below the critical densities of the higher frequency 
lines for high dipole moment molecules. For these regions 
statistical equilibrium (SE) calculations, also referred to 
as non-LTE or NLTE, need to be performed. 

The populations of the levels are calculated by solving 
the equation: 

Y^[njAji + {rijBij - mBi^Jji] 

3>l 

- Yl\ niAl i + ^ UlBl i ~ n i B iOJij] (2) 

3<l 

+ Yy>>jCji - niCij] = 0, 
j 



where Jji is the integrated mean intensity, Aij and Bij 
the Einstein coefficients and Cij the collisional rates. The 
Einstein A coefficients and collisional rate coefficients for 
CO, HCO + and HCN are the same as those in Jansen et al. 
(1994, Table 4). The calculation of the level populations is 
an iterative process since the integrated mean intensity is 
directly related to the levels ni, which in turn affects the 
mean intensity. 

Test calculations have been performed for three cases. 
The first is LTE, where the populations are given by 
the Boltzmann equation. In this assumption the pop- 
ulations are dominated by collisions and therefore de- 
pend only on the local temperature. The second is 
Statistical Equilibrium without stimulated radiative ef- 
fects (SE[J„=0]), in which the populations of the lev- 
els are no longer assumed to be dominated by collisions 
and are calculated explicitly. In this case, Equation 
is solved under the assumption that J u —0 for all radia- 
tive transitions. The third method is the full Statistical 
Equilibrium (SE) solution using a Monte Carlo code 
( Hugcrhcijdc fc van dcr Tak 200C| ) to calculate the mean 
intensity J u for each radiative transition iteratively, tak- 
ing the line emission and absorption throughout the 2D 
disk into account. 

In Figure [|, the relative populations for the levels J— 3 
and 6 of CO and J—l and 4 of HCO + arc plotted for the 
LTE, SE(7„=0) and SE for the D'Alessio et al. (1999) 
model. This model is chosen because it has a smooth tem- 
perature gradient but does show a temperature inversion 
in the ^-direction. The adopted CO and HCO + abun- 
dances are 10 and 5T0 -9 respectively, and the turbulent 
line width is assumed to be 0.2 km s _1 . Figure ^ shows 
that the differences between the three methods are small 
in the midplane but that they become significant in the 
lower density upper layers. For the J=l level of HCO + , the 
influence of the Cosmic Microwave Background 2.7 K ra- 
diation is apparent since its relative population continues 
to rise toward the outside in the SE calculation compared 
to the SE(J„=0) calculation. 

Even though only the full 2D SE calculation describes 
the populations accurately, its calculation is an enormous 
computational task due to the large column densities and 
steep gradients in density and temperature coupled with 
a narrow velocity profile. In these cases, the convergence 
criteria of a numerical code become very important. The 
SE( J u =0) calculation provides better agreement with the 
SE calculation compared to LTE, especially at larger dis- 
tances from the star where most of the observed radiation 
originates. Therefore the SE(J„=0) method is adopted in 
the ray-tracing calculations. 

The abundances and depletion of various molecules are 
taken into account in two different ways: 

1. A constant model, in which the fractional abundances 
are the same throughout the disk. In the standard 
model, abundances close to those found in dark clouds 
are chosen. These values can be subsequently lowered 
by a constant factor Dq with respect to the standard 
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Fig. 5. The two ways in which the depletion of molecules com- 
pared to the interstellar value is taken into account in the 
computations. On the left a constant depletion Da is shown, 
whereas on the right a step-function representing the freeze-out 
of molecules Dj below a critical temperature is illustrated 



values. Note that the chemical interpretation of Dc 
may vary: it can be due to freeze-out, photodissoci- 
ation, different chemistry at high or low densities, or 
any combination of these. A sketch of the depletion 
models compared to interstellar abundances is given 
in Figure ||. 

A jump model, in which the abundances are constant 
except for regions where the temperature is lower than 
a critical temperature. Below this temperature the 
abundance is assumed to drop by a certain factor due 
to freeze-out of the molecule onto cold dust grains. In 
this paper, a temperat ure of 22 K (based on CO freez- 
ing onto a CO surface; Bandford fc Allamandola 1993 ) 
is assumed (see Figure ^|). HCO + is assumed to fol- 
low the CO abundance profile and will thus deplete 
at the same temperature. For HCN a critical tempera- 
ture of 80 K is assumed. A more detailed and realistic 
description of the depletion and abundance variations 
has been given by Aikawa & Herbst (1999) and Willacy 
et al. (1998), but these results are too specific for the 
exploratory purposes of this paper. 



4.4. Approximate line intensities 
4.4.1. Vertical radiative transfer 

To constrain the depletions and thereby bracket the com- 
putational domain of the problem, the intensity of each 
line is calculated by solving the radiative transfer equa- 
tion 



— = a v { — 

as a v 



(3) 



in the vertical direction, where l v is the intensity, s the 
path length along a ray normal to the disk, j v the emission 
function, and a v the absorption function, which is the 
inverse of the mean free path. The ratio of the emission 
and absorption functions is known as the source function. 
The transfer equation is solved in an iterative ray-tracing 
procedure using small As steps from sq = — oo to sjv = 
+oo (the observer). 



The equation to be solved thus becomes 

-feline (Si) = Iv{Si-l)e~'' 



Ju {si)<t>{vy '-As, 



(4) 



where <j){v) is the line emission profile function in fre- 
quency space and I v {sq) = 0. 

The optical depth t v is the sum of the attenuation by 
dust and gas along the step-length and is equal to 



(a c + ai) As 



where 



Oil 



Rgd 
AulC 3 



9u 
xi — 

9i 



8ttv 3 

and the emission is given by 



x u n gas A m </>(i/) 



jv = -7-A u ix u n sas X m <j)(v). 



(5) 

(6) 
(7) 

(8) 



In these equations, m gas is the mean molecular weight of 
a gas particle, R g( j the gas to dust ratio, A u i the Einstein 
coefficient, v the frequency of the transition, c the speed 
of light, X m the abundance of the molecule relative to 
the gas density n gas , which is taken to be equal to the 
density of H2, gi the statistical weight of level i and x\ 
the population of level i. A mean molecular weight of 2.2 
(proton mass per particle) and R g d = 100 are adop ted . 
The level populations are solved using SE (J v =0) (§4.3) 
for the reasons stated above. 

The continuum mass absorption coefficient n u is taken 
from Osscnkopf fc Hcnning (1994) and extended to wave- 
lengths longer than 1.3 mm as 



k c (1.3mm) 



2.31 x 10 11 Hz 



(9) 



where K c (1.3mm) depends on the specific coagulation 
model chosen and (3 is taken to be 1. For this problem 
the k c values were taken from the model with coagula- 
tion at a density of 10 s cm -3 covered by a thin ice-layer 
(k c (1.3 mm)=1.112 cm 2 g _1 ). In practice, however, the 
dust absorption is negligible compared to the line absorp- 
tion for the molecular transitions in the wavelength range 
of interest. 

4.4.2. Calculation of intensity 

For comparison with observations and model results, all 
values are referred to the size of the disk model. The ob- 
servations are thus scaled as follows: 



<Jisk,obs 



n 



!/,tel 



9. 



disk 



^tel °„,disk 

2tt J^ max RdR 1 



AIT 



(rf(pc)) 2 



(10) 



(11) 
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where the ratio T m b = T^/rj^ is given in Table [l], ^„tel 
is the beam size at the frequency of the line in arcsec 2 , 
d the distance in parsec, AU the astronomical unit, and 
^disk ^ ne s i ze °f the disk in arcsec 2 . 

The models are calculated using the ray-tracing 
Equation (||) which are weighted according to its emitting 
surface area 



RdR 



(12) 



Here, / is the mean intensity for the entire disk at the 
surface, with Nx the number of cells in the R— direction. 
The disk is assumed to be seen face-on. This is a very 
good approximation for TW Hya, whereas it should still 
be reasonable for LkCa 15 even though the disk is seen 
at an inclination of 60°. The intensities derived by this 
method will be used only in the analysis of line ratios, 
which are less sensitive to inclination effects than absolute 
values. 

The radius of the disk is taken to be 200 AU in all 
models, or 400 AU diameter. The size of the LkCa 15 disk 
(d=140 pc) suggested by the OVRO 13 CO maps of Qi 
(2000) is slightly larger (420 AUx530 AU). For TW Hya, 
no millimeter interferometer observations are available, 
but mid-infrared and VLA 7 millimeter images suggest a 
disk size of -100 AU flWilncr ct al. 200C| ). Scattered light 
images observed with the Hubble Space Telescope suggest 
an outer radius of at least 200 AU, however (Weinberger 
et al. 1999, Krist et al. 2000). We therefore adopt a similar 
disk size in AU as for LkCa 15, but with d=57 pc. 

4.5. Line intensities using full radiative transfer 

For the calculation of the populations in SE using full 
radiative transfer, the Monte Carlo code developed by 
Hogerheijde & van der Tak (2000) is used. In this code, 
Equation (^) is solved in an iterative fashion, where all 
photons start at the outer boundary with an intensity 
given by the 2.728 K Cosmic Background radiation. In this 
calculation, the inferred abundances from the SE(J„=0) 
method are adopted. The calculated populations at each 
position in the disk are used to compute the complete line 
profiles of selected molecules using a program which cal- 
culates the sky brightness distribution. The profiles are 
calculated by constructing a plane through the origin of 
the disk perpendicular to the line of sight, with a spatial 
resolution small enough to sample the physical and ve- 
locity distributions. Both the spatial resolution and the 
velocity resolution can be specified. A ray-tracing calcu- 
lation is performed through this plane from — oo to +oo, 
keeping track of the intensity in the different velocity bins. 



5. Results 

The models are calculated initially using standard dark 
cloud abundances of CO of M0~ 4 , HCO+ of 5-1Q- 9 , 



and HCN of 5T0~ 9 relative to H2. For the isotope ra- 
tios the following values are used throughout: 12 C/ 13 C=60 
and 16 O/ 18 O=500. These models are referred to as 
Dc=Dj=l. Subsequently, the values of Dc and Dj are 
varied (see §4.3). The line intensities have been calculated 
assuming a micro-turbulence of 0.2 km s _1 . For TW Hya, 
this results in calculated line widths of ~ 0.8 km s _1 , in 
good agreement with observations. 

5.1. t = 1 surfaces 

Significant insight into the observational results can be 
obtained by investigating the regions of the disk where 
the different lines become optically thick. At each radius 
the effective emission region for each line is calculated us- 
ing the SE(J„=0) method by integrating from the top 
layer down until r=l in line + continuum is reached. 
Although the r=l level is chosen arbitrarily and radia- 
tion from deeper in the disk may still escape, it provides 
a useful measure of the volume of the emitting region for 
each molecular transition. This calculation is performed 
only for a face-on disk for simplicity and is thus only ap- 
plicable for the TW Hya case. It does, however, give an 
indication of the parts of the disk from which the molecu- 
lar emission arises in more general cases. For the model by 
D'Alessio et al. (1999), a contour-plot of the r = 1 surfaces 
of the observed CO and HCO + lines is given in Figure |[ 
where the former are overplotted on the temperature dis- 
tribution and the latter on the density distribution. The 
line emission is dominated by densities and temperatures 
above the r = 1 contour, which can then be compared 
to the values derived from the constant temperature and 
density models given in §^|. 

It is seen that, for standard abundances, the 12 CO lines 
become optically thick in the upper, warm layer of the 
disk where T > 40 K. On the other hand, the 13 CO lines 
probe into the colder regions around 20-30 K. Thus, the 
12 CO excitation temperature should be higher than that of 
13 CO, which must be taken into account in the analysis of 
isotopomeric line ratios. Similarly, the higher frequency 3- 

2 and 6-5 lines generally have higher optical depths than 
the 1-0 lines, and thus probe better the warmer upper 
layer. Even C ls O is not fully optically thin, but has r ss 
1-2. The low temperature of 20-30 K probed by 13 CO 
is consistent with the simple analysis of the data in §3. 

For the standard HCO + abundance, the 1-0 to 4- 

3 lines are optically thick in the outer layers, whereas 
the H 13 CO + lines are close to optically thin throughout 
the disk. Thus, the HCO + lines probe densities of order 
10 6 — 10 7 cm~ 3 , below the critical density of the 4-3 tran- 
sition. For H 13 CO + , the populations will be closer to LTE 
because its emission arises primarily from regions with 
densities of 10 7 — 10 s cm -3 . If the HCO + abundance is 
depleted by a constant factor Dc=10, the HCO + lines be- 
come optically thin in the outer regions of the disk and 
now trace regions with densities of 10 7 — 10 8 cm -3 . The 
HCN 1-0 to 4-3 lines show a similar behavior to HCO + . 
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Fig. 6. The r = 1 surlaces for the observed CO and HCO + isotopomer lines integrated from the top, overplotted on the 
temperature (top) and density distribution (bottom) in the model by D'Alessio et al. (1999). The dotted contours are iso- 
temperature (in K) or iso-density (in cm -3 ) contour lines. The results are shown for both the standard abundances (left) and 
depleted by a constant factor D c = 10 (right). The labels indicate: A: 12 CO 6-5; B: 12 CO 3-2; C: 12 CO 2-1; D: 13 CO 3-2; E: 
13 CO 1-0; F: C ls O 3-2; G: C ls O 2-1; H: HCO+ 4-3; I: HCO+ 1-0; J: H 13 CO+ 4-3; K: H 13 CO+ 1-0. 



The densities of 10 6 — 10 s cm -3 derived from the observed 
HCO + and HCN lines in §3 are consistent with this anal- 
ysis for modest depletions of both species. 



5.2. Integrated line ratios: range of depletions 

In this section, the relative line intensities obtained in the 
SE( J„=0) method are used to constrain the abundances of 
the molecules and the level of depletion. Since lines of dif- 
ferent isotopomers arise from different regions, their line 
ratios will depend on the local depletion values. By cal- 
culating models for a range of depletions, the abundances 
can be derived by varying both the overall de plet ion Dc 
and the jump depletion Dj as described in §4.2. In the 



comparison of the line ratios, the difference in beam dilu- 
tion for the two lines must be taken into account (§ 4.4.2|) . 



Isotope ratios are more sensitive to both Dc and Dj 
as compared to ratios of different species. For instance the 
12 CO emission remains optically thick up to large values 
of Dc and therefore does not probe the region below 22 K, 
whereas 13 CO becomes sensitive to Dj for modest values 
of D c - In Figure @ the CO/ 13 CO 3-2 and 13 CO/C 18 
3-2 line intensity ratios are plotted as functions of Dj 
and Dc- The observed values for LkCa 15 are plotted 
with dashed lines, and indicate a range Dc G [3, 40] when 
both plots are combined. The ratio for TW Hya (dot- 
dashed line) indicates a larger depletion with Dc > 100 
and Dj > 10. Combining similar plots for all species and 
lines, the resulting values of Dc and Dj are shown in 
Table || for the three models of interest for both sources. 
The inferred ranges for the disk models are large and it is 
difficult to give accurate values for molecules for which few 
lines have been measured. The abundances of all molecules 
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D'Alessio 


et al. (1999) 


Chiang & 


Goldreich (1997) 


Bell (1999) 






LkCa 15 


TW Hya 


LkCa 15 


TW Hya 


LkCa 15 


TW Hya 


CO 


D c 


[3,15] 


>30 


[3,30] 


> 100 


[1,500] 


[10,1000] 




D., 


[3,30] 


>1 


[1,15] 


>1 


>1 


>1 


HCO+ 


D c 


[3,80] 


> 80 


[10,100] 


> 100 


[1,80] 


[2,1000] 




Dj 


>1 


>1 


>10 


>1 


[1,100] 


>1 


HCN 


D c 


[1,400] 


[4,600] 


[10,200] 


[10,800] 


[1,500] 


[4,500] 




Dj 


>1 


>1 


>1 


>1 


>1 


>1 



The numbers in square brackets indicate the range of inferred depletions 
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Fig. 7. The CO/ 13 CO 3-2 (top) and 13 CO/C 18 3-2 
(bottom) line intensity ratios as functions of the jump 
depletion Dj and the overall depletion Dc within the 
D'Alessio et al. (1999) model. The observed ratios for 
LkCa 15 CO 3-2/ 13 CO 3-2 (dashed line) and TW Hya 
(dash-dotted line) are shown in the figures. 



5.3. Line profiles and intensities 

The line profiles are calculated using the full 2D radia tive 
transfer code for the range of depletions derived in § |5.2| . 
The depletions are further constrained by the absolute 
intensities. Specifically, for LkCa 15 Z?c=5 and 10 with 
Dj=10 is taken, and for TW Hya the same Dj=10 was 
assumed but with _Dc"=100 and 200. As a reference, an 
extra run was performed for LkCa 15 with no depletions. 
A general turbulent width of 0.2 km s _1 is assumed and 
the only structured velocity distribution is taken to be the 
Kcplcrian rotation of the disk. This velocity component is 
important for the comparison with observations of sources 
at non-zero inclination. For these calculations, an inclina- 
tion of 60° for LkCa 15 and 0° for TW Hya is used. The 
results are convolved with the appropriate beam as given 
in Table |l|. 

A model with no depletion was also run for LkCa 15 
with an inclination of to check the effect of inclina- 
tion. Although the total integrated line intensities changed 
significantly, their ratios changed only up to 7 %. This jus- 
tifies the approximate radiative transfer approach used in 
§5.2 for a first estimate of the depletions. 

The resulting integrated intensities are presented in 
Table 3 for the high- J rotational lines and in Table 4 for 
the lower- J transitions. For six high- J rotational lines, the 
observed profiles are plotted in Figure || with the three 
calculated model emission profiles superposed. In the left- 
hand figures, three lines are shown for LkCa 15 whose clear 
double peaks are due to the Keplerian rotation in the disk. 
On the right, the single peaks for a face-on disk such as 
that around TW Hya are seen. The optically thick lines 
from the latter source show that the disk can be fitted 
with a turbulent velocity of 0.2 km s . 



5.3.1. Depletions 



in the TW Hya disk seem to be depleted by a large factor. 
In general Dj w 10 is taken as a best fit for both sources. 
In cases where no constraints are available, Dj = 10 has 
been assumed. 



The absolute intensities in Tables 3 and 4 indicate that 
refinements of the inferred depletions are required, since 
different molecules favor different amounts of depletion. 
Note that the intensities computed for the cold Bell model 
are always smaller compared to the other two models, for 
both the high and low rotational lines. The reason for this 
is twofold. First, for a cold isothermal disk structure the 
level populations of any molecule at densities above the 
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Table 3. Integrated intensities for the higher rotational lines for the three disk models 
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u. iy 




U.UOl 




9/1 








LkCa 15 d 




0.53 




1.39 


0.39 


<0.14 


0.25 








0.26 




<0.13 




Model 


Dc 
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CO 6- 
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CO 3-2 


"CO 3-2 


HCN 4-3 
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2 
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3 


HCO + 4- 
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A 
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10 
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0.29 


0.36 
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A 
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10 
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C 


200 


10 


0.57 




1.73 


0.32 


0.43 


0.21 




0.025 




0.55 




0.040 




B 


10 


10 


0.38 




1.41 


0.15 


0.54 


0.48 




0.045 




0.62 




0.06 




TW Hya d 




<3.2 




1.98 


0.24 


0.49 


0.45 




<0.04 




0.49 




0.07 





a D'Alessio et al. (1999) model 

b Bell (1999) model 

c Chiang & Goldreich (1997) model 

d The observed values have an estimated uncertainty of 20 %; all values refer to the original beam size of the observations (see 
Table |l|) 



Table 4. Integrated intensities for the lower rotational lines for the three disk models 



Model Dc 


Dj 


13 CO 1-0 


C 18 1-0 


HCN 1-0 


H 13 CN 1-0 


HCO+ 1-0 


H 13 CO+ 1-0 




1 


1 


6.18 


1.90 


4.47 


0.85 


4.84 


0.20 


B 6 


1 


1 


4.04 


1.57 


2.58 


0.76 


2.68 


0.17 


C c 


1 


1 


8.80 


2.88 


5.72 


0.96 


5.94 


0.25 


A 


5 


10 


3.60 


0.53 


2.70 


0.11 


3.11 


3.3E-2 


B 


5 


10 


2.04 


0.32 


1.39 


5.5E-2 


1.56 


1.8E-2 


C 


5 


10 


3.56 


0.62 


2.52 


0.14 


3.06 


4.0E-2 


A 


10 


10 


2.31 


0.31 


1.75 


6.0E-2 


2.11 


1.8E-2 


B 


10 


10 


1.55 


0.18 


1.11 


2.8E-2 


1.31 


9.3E-3 


C 


10 


10 


2.40 


0.38 


1.70 


7.6E-2 


2.10 


2.3E-2 




LkCa 15" 




7.43 


0.70 


3.04 


1.20 


4.19 


7.E-2 



a D'Alessio et al. (1999) model 

b Bell (1999) model 

c Chiang & Goldreich (1997) model 

d The observed values, taken from Qi (2000), have an estimated uncertainty of 20 %; all values refer to the original beam size 
of the observations (see Table 0) 



critical density are the same at each position. This means 
that the optical depth becomes directly proportional to 
the column of gas. For a model with an increasing temper- 
ature, the optical depth will not be directly proportional 
to the column but smaller. A model with a step function 
in its temperature will give results in between these two 
cases. Second, the colder disk will have slightly narrower, 
more optically thick lines due to the lower thermal motions 
in the gas, thereby trapping radiation more effectively. 

Both effects are visible in the 13 CO and CO lines (Fig. 
||). The two peaks in the 13 CO intensity for the inclined 
LkCa 15 disk are reduced significantly compared to the 



emission at line center in the Bell model. In the CO 3-2 
line for the face-on TW Hya disk, the emission predicted 
by the Bell model is extremely optically thick, shown by 
the flat-topped emission profile. Also, the total lincwidth 
is somewhat smaller compared to the other two disk mod- 
els due to the low temperatures. Thus, the observed line 
profiles argue for a rising temperature structure in the ver- 
tical direction to prevent the high optical depths found in 
the cold isothermal model. 

To counteract the low intensities found in the Bell 
model for the TW Hya disk, additional calculations were 
performed for less severe depletions (Table 3: Dj=10, 
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0.20 




-6 -2 2 6 -2 -1012 




Fig. 8. Observed line profiles for LkCa 15 (left) and 
TW Hya (right) with the three models superposed (solid: 
D'Alessio et al.; dashed: Chiang & Goldreich; dash-dotted: 
Bell). The different transitions are indicated with the 
adopted depletions compared to standard molecular cloud 
values. 



Dc = 10). The integrated intensities increase to just above 
the observed values in this case; however, the lines are ex- 
tremely optically thick and show nearly square emission 
profiles, which is not observed for the CO 3-2 and HCO + 
4-3 lines. Thus, the line profiles speak against small deple- 
tions. In the two warm disk models, there is no significant 
difference between the two assumed depletions and only a 
slight preference can be given to Dq—WO. This is largely 
based on the HCO+ /H 13 CO+ ratio, which is ill fitted by a 
depletion of 200 and only moderately well for Dc = 100. 

For LkCa 15, CO is best fitted with little depletion: all 
three models indicate Dc close to 1 for the lower rotational 
lines. The upper limit on the C 18 3-2 line and the C 18 
1-0 emission would favor some depletion, which, in the case 
of C 18 0, can be explained by enhanced photodissociation 
in the upper layers due to the lack of self shielding. The CO 
3-2 intensity is too low in all three models, which may be a 
sign of extended emission beyond 200 AU since this line is 
optically thick. The HCO + and HCN lines are best fitted 
by a moderate depletion of Dc=5 and Dj=10. For HCN, 
this could again be a sign of a lack of shielding against 
photodissociation compared to CO. The observed H 13 CN 
1-0 is slightly too high for all three models. Together, the 
HCN and H 13 CN data indicate that the HCN abundance 
needs to be lowered primarily in the surface layer to both 
bring the main isotope HCN emission down but keep a 
high H 13 CN intensity. 



5.3.2. Probing the temperature and density structure 

The calculated line ratios which are sensitive to the tem- 
perature and density distribution are summarized in Table 
[|. The temperature of the upper layer in the LkCa 15 
models, as probed by the CO 6-5/3-2 ratio, fits within 
the errors to all three models, confirming the rather cold 
upper layer of this disk. However, as explained in the pre- 
vious section, the absolute intensities are too low in the 
cold Bell model. To reproduce the observed CO 6-5 inten- 
sity, the CO abundance would have to be increased well 
above the cosmic carbon abundance in the cold isothermal 
model. 

For the TW Hya disk, the modeled CO 4-3/3-2 ratios 
are on the low side, even in the D'Alessio et al. and Chiang 
& Goldreich models, indicating that the temperature in 
the surface layers would need to be higher. However, the 
calibration uncertainties in the older Kastner et al. (1997) 
data make this conclusion less firm. Further observations 
of high- J CO lines are needed to constrain the tempera- 
ture structure of this disk. 

The density is probed by the different HCO + and HCN 
ratios. All three models are consistent with the observed 
4-3/1-0 ratio for the main isotopes, indicating that the 
density in the layer above the midplane is in the correct 
range. The upper limits on the H 13 CN and H 13 CO + pre- 
vent any conclusions for the midplane. As Table ^ shows, 
the models make different predictions for these optically 
thin species and future observations may be used to dis- 
tinguish them. 

Overall, the absolute line intensities and ratios are con- 
sistent with the models of D'Alessio et al. and Chiang & 
Goldreich for reasonable values of the depletions. The cur- 
rent data cannot distinguish between these two flared disk 
models. There is some evidence, however, both from the 
line ratios and from the line profiles that the surface of the 
disks needs to be warmer than that of a shielded isother- 
mal outer disk such as computed by Bell (1999). 

6. Discussion 

The high depletions derived for TW Hya are in agreement 
with the conclusions by Kastner et al. (1997). The deple- 
tion can be due to two reasons: the first is a change in 
the gas- to dust-mass ratio to a value lower than 100 due 
to removal of gas. The second possibility is a large deple- 
tion of CO and other molecules, but not H2. The former 
possibility could be partly tested by searches for the pure 
rotational H2 lines (Thi et al. 2001). Regarding the sec- 
ond option, the present analysis indicates that the deple- 
tion cannot simply be due to the freezing out of molecules 
resulting in a large value of Dj: the molecules also need 
to be depleted in the warmer upper layers probed by Dc- 
This latter conclusion could be consistent with the fact 
that TW Hya is a very active UV and X-ray producing 
star (Kastner et al. 1999), capable of destroying CO due to 
dissociation or ionization. This can be tested by measuring 
the CO ionization or dissociation products, such as C + , C 
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Table 5. Line ratios obtained for the three disk models 



Line ratio 


Observed 


D'Alessio et al. 


Bell 




Chiang & 


I Goldreich 




LkCa 15 c 


D c , Dj 


Ratio 


D c , Dj 


Ratio 


D C , Dj 


Ratio 


CO §5§ 




5,10 a 


0.66 


1,1° 


0.36 


10,10 


0.55 


13 co f=§ 


0.05lH 3 2 


1,1° 


0.06 


10,10 a 


0.05 


1,1° 


0.05 


HCO+ |5§ 


0.06lH 3 2 


1,1 


0.11 


1,1° 


0.06 


1,1 


0.08 


H 13 CO+ 4=§ 


< 1.9 


10,10 6 


1.44 


5,10 


1.00 


io.io 6 


1.78 


HCN |5§ 


0-081^ 


5,10 a 


0.09 


5,10 a 


0.07 


l,l a 


0.07 


TW Hya c 


CO §E| 


< 1.62 


100,10 


0.48 


100/10,10 d 


0.12/0.27 


100,10 


0.41 


CO|5| 


o cq+1.26 
i - JJ -0.84 


200,10 


1.28 


200/10,10 d 


1.01/1.13 


100,10 


1.29 


HCN |=| 


1 09 +0M 


100,10 


1.64 


100/10,10 d 


1.38/1.13 


100,10 


1.92 



a Ratios for all three combinations of Dc and Dj fall within the error 
6 Ratios for D c = 5, Dj = 10 and D c = 10, Dj = 10 fall within the error 

c The observed values and all ratios refer to the original beam sizes in which the lines were observed (see Table |l|); thus, the 
beam size differs between species and lines. The error bars correspond to a 20 % uncertainty. 
d Ratios for high and low values of Dc are shown 



and CO + . TW Hya seems well described by a flared disk 
model where the ultraviolet radiation is capable of heating 
the upper layer (see §5.3). 

The modeling suggests freezing out with a common 
value of Dj of at least 10. The depletion of molecules 
onto grains can be an important chemical sink in these 
disks since most of the mass is cold. This makes mass 
determinations of disks using CO or any of its isotopomers 
highly uncertain. 

New models have recently been fitted to the SEDs of 
LkCa 15 and TW Hya by Chiang et al. (2001) and Chiang 
(2000). These models use as one of the main parameters 
the dust settling toward the midplane. For LkCa 15, the 
SED modeling indicates that the dust should have settled 
within the disk scale height H to explain the observations, 
with a high dust temperature in that region (Tdust =49 K 
at 100 AU). In these high density regions, the gas and dust 
temperature should be coupled; however at heights above 
the scale height H the lack of grains will reduce the heat- 
ing of the gas due to the photo-electric effect, although 
some small grains may still be present. Together with en- 
hanced cooling due to [O I] and [C II], this may cause a 
second temperature inversion with a cool upper layer free 
of dust. The lack of grains in the upper layer would be con- 
sistent with the non-detection of LkCa 15 by HST optical 
images (K. Stapelfeldt, private communication). This sug- 
gests a lack of scattering which should have been readily 
seen for a flared disk at an inclination of 60° with grains 
well mixed with the gas and an albedo of 0.5. In addition, 
the models used were calculated for a stellar temperature 
of 4000 K, whereas LkCa 15 has a higher effective tem- 
perature (T e g:=4400 K) which would result in higher disk 
dust-temperatures. The relatively low gas temperature in 
the disk, indicated by the line observations, strengthen 
the conclusions derived from the scattering and SED ob- 



servations that dust settling has taken place in LkCa 15. 
Self-consistent models of the gas temperature and abun- 
dances of LkCa 15 are needed, taking dust-settling into 
account. 

As noted above, it is not yet possible to distinguish 
between the D'Alessio et al. (1999) and the Chiang & 
Goldreich (1997) models with the current observations. 
The data lack spatial resolution and have insufficient sen- 
sitivity to observe the optically thin isotopic lines. In ad- 
dition to higher spatial resolution and sensitivity, better 
calibration of the data is needed, all of which will be pro- 
vided by the Atacama Large Millimeter Array. 

7. Conclusions 

The main conclusions from this work are: 

— High-frequency molecular lines with high critical den- 
sities and excitation temperatures are detected from 
circumstellar disks. These observations can be used to 
test the temperature and density structure of different 
disk models in the literature. 

— The t — 1 surfaces of the various lines indicate that the 
observed emission of the main isotopes originates from 
the (warm) intermediate layer of the disk, whereas the 
emission from the 13 C isotopes may also probe the mid- 
plane if the molecule is not frozen out completely. 

— Most molecules are depleted by a large factor (> 100) 
for TW Hya and a smaller factor (surface w 1-5, mid- 
plane fa 1 — 50) for LkCa 15. Freeze-out onto grains 
at T < 22 K is indicated by the observations, but the 
molecules are also depleted in the upper, warmer lay- 
ers, likely due to photodissociation. Gaseous species 
like CO and its isotopomers should therefore not be 
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used as mass tracers due to their uncertain abun- 
dances. 

— A model with a cold isothermal temperature distribu- 
tion will have high optical depths in the lines, thereby 
reducing the integrated line emission and producing 
flat-topped profiles. 

— The TW Hya disk has a warm surface layer (>40 K) 
while the LkCa 15 surface layer is cooler (~20-40 K). 
This conclusion depends sensitively on the calibration 
accuracy of the high- J CO lines. 

— The inferred warm dust from the SED combined with 
the cooler gas detected here may be consistent with 
settling of dust in the LkCa 15 disk. 

— The density profiles in the three models are consistent 
with the observed line ratios, except that the temper- 
ature in the upper layer of a non-irradiated disk such 
as in the Bell (1999) model is too low for some sources. 
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